##########################################################################################################
## Payson (2021) When Cities Lobby Replication Code
## Chapter Four: Figures (4.1 - 4.7) 
## Created with R version 4.0.3 (2020-10-10) -- "Bunny-Wunnies Freak Out"
##########################################################################################################


## Load packages
library(maps); library(ggplot2); library(scales)
library(RColorBrewer); library(tidyverse); library(ggpubr)


## Set working directory
setwd()



##########################################################################################################
## Figure 4.1 Proportion of Cities Lobbying by State
##########################################################################################################


dta <- read.csv("city_lobby.csv")
states_map <- map_data("state")


## Subset to balanced data
dta <- dta[dta$year > 2005,]


## Get total # of cities and % cities ever lobbying by states
plot.dta <- dta %>%
  group_by(state) %>%
  summarize(ever.lobby = sum(ever.lobby) / 10,
            num.city = length(state) / 10,
            x = mean(long),
            y = mean(lat))


plot.dta$pct.lobby = (plot.dta$ever.lobby / plot.dta$num.city) * 100
plot.dta <- subset(plot.dta, state != "Hawaii" & state != "Alaska")
plot.dta$state <- tolower(plot.dta$state)

## Create state lobbying categories
plot.dta$pct.lobby1 <- ifelse(plot.dta$pct.lobby < 15, 0,
                               ifelse(plot.dta$pct.lobby < 30, 15,
                                      ifelse(plot.dta$pct.lobby < 50, 30,
                                             ifelse(plot.dta$pct.lobby >= 50, 50, NA))))

add.white <- plot.dta[plot.dta$state == "california" | plot.dta$state == "arizona", ]


(p <- ggplot(plot.dta, aes(x=x, y=y, map_id = state)) +
    geom_map(map = states_map, colour = "black", aes(fill = factor(pct.lobby1))) +
    geom_map(data = add.white, map = states_map, colour = "white", fill = "black") +
    scale_fill_manual(name = "% Cities Lobbying",
                      values = c("white", "gray70", "gray40", "black"), 
                      labels = c("0-15%", "15-30%", "30-50%", "Over 50%")) +
    theme_minimal(base_family = "serif", base_size = 10) +  
    xlim(-124, -69) + ylim(25, 49) +
    theme(plot.background = element_blank(),
          legend.position = "top",
          legend.key.width = unit(1, "cm"),
          axis.title.y = element_blank(),
          axis.text.y = element_blank(),
          axis.ticks.y = element_blank(),
          axis.title.x = element_blank(),
          axis.text.x = element_blank(),
          axis.ticks.x = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank()))



##########################################################################################################
## Figure 4.2 Geographic Correlates of City Lobbying Across States
##########################################################################################################

dta <- read.csv("state_data.csv")

dta$termlim <- factor(dta$termlim)
dta$year <- factor(dta$year)
dta$limits <- factor(dta$limits)


## Specify model
m <- lm(lobby.pct ~ limits + city.density + frag + termlim*leg.prof + year, data = dta)


## Get predicted values
pred <- predict(m)

plot.dta <- data.frame(pred = pred,
                       state = dta$abbr,
                       city.density = m$model$city.density,
                       limits = m$model$limits,
                       frag = m$model$frag,
                       termlim = m$model$termlim,
                       leg.prof = m$model$leg.prof,
                       year = m$model$year)


## Predicted probability by average city population
(p1 <- ggplot(plot.dta, aes(x = city.density, y = pred)) +
    stat_summary_bin(geom = "point", bins = 50) +
    geom_smooth(method = lm, formula = y ~ poly(x, 2), color = "black", se = FALSE) +
    labs(y = "% Cities Lobbying", x = "Average Residents Per City (Log)")  +
    theme_bw(base_family = "serif", base_size = 11) +
    theme(plot.background = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank())) 


## Predicted probability by city fragmentation
(p2 <- ggplot(plot.dta, aes(x = frag, y = pred)) +
    xlim(0, 8) +
    stat_summary_bin(geom = "point", bins = 50) +
    geom_smooth(method = lm, formula = y ~ poly(x, 1), color = "black", se = FALSE) +
    labs(y = "% Cities Lobbying", x = "Cities Per Thousand Square Miles of Land")  +
    theme_bw(base_family = "serif", base_size = 11) +
    theme(plot.background = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank())) 


(fig <- ggarrange(p1, p2, ncol = 2, nrow = 1))


##########################################################################################################
## Figure 4.3 Local Tax Limits and City Lobbying Across States
##########################################################################################################

dta <- read.csv("state_data.csv")

dta$termlim <- factor(dta$termlim)
dta$year <- factor(dta$year)
dta$limits <- factor(dta$limits)


## Specify model
m <- lm(lobby.pct ~ limits + city.density + frag + termlim*leg.prof + year, data = dta)


## Get predicted values
pred <- predict(m)

plot.dta <- data.frame(pred = pred,
                       state = dta$abbr,
                       city.density = m$model$city.density,
                       limits = m$model$limits,
                       frag = m$model$frag,
                       termlim = m$model$termlim,
                       leg.prof = m$model$leg.prof,
                       year = m$model$year)


(p <- ggplot(plot.dta, aes(x = limits, y = pred)) +
   geom_boxplot(outlier.shape = NA, coef = 1.25, notch = T) +
   labs(y = "% Cities Lobbying", x = "Number of Property Tax Limitations")  +
   theme_bw(base_family = "serif", base_size = 12) +
   theme(plot.background = element_blank(),
         legend.position = "top",
         panel.grid.major = element_blank(),
         panel.grid.minor = element_blank())) 



##########################################################################################################
## Figure 4.4 State Legislative Professionalism, Term Limits, and City Lobbying
##########################################################################################################

dta <- read.csv("state_data.csv")

dta$termlim <- factor(dta$termlim)
dta$year <- factor(dta$year)
dta$limits <- factor(dta$limits)


## Specify model
m <- lm(lobby.pct ~ limits + city.density + frag + termlim*leg.prof + year, data = dta)


labels <- plot.dta %>%
  group_by(state) %>%
  summarise(pred = mean(pred),
            leg.prof = mean(leg.prof))

labels <- subset(labels, state=="CA" | state=="MI" | state == "OH" | state == "AZ" | state == "NV" |
                   state == "CO" | state == "FL" | state == "TX" | state == "ND" | state == "NH" |
                   state == "DE" | state == "NJ" | state == "IL" | state == "PA" | state == "MA" |
                   state == "WI" | state == "NY" | state == "NC" | state == "TN" | state == "TN")


(p <- ggplot(plot.dta, aes(x = leg.prof, y = pred, linetype = termlim)) +

    stat_summary_bin(geom = "point", aes(color = termlim), bins = 100) +
    geom_smooth(method = lm, formula = y ~ poly(x, 2), aes(color = termlim), se = FALSE) +
    ylim(0, 40) +
    scale_linetype_manual(values = c("solid", "dashed")) +
    scale_colour_grey(start = 0.6, end = .2,
                      name = "", labels = c("No Term Limits", "Term Limits")) +
    labs(y = "% Cities Lobbying", x = "Legislative Professionalism")  +
    theme_bw(base_family = "serif", base_size = 12) +
    guides(linetype=FALSE) +
    theme(plot.background = element_blank(),
          legend.position = "top",
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank()) +
    annotate("text", x = labels$leg.prof, y = labels$pred + 2, label = labels$state, size = 4))



##########################################################################################################
## Figure 4.5 Political Correlates of City Lobbying
##########################################################################################################

dta <- read.csv("state_data.csv")

dta$gov.partyR <- ifelse(dta$gov.party == "R", 1, 0)
dta$leg.controlR <- ifelse(dta$leg.control == "R", 1, 0)
dta$house.rep.pct <- dta$house.rep.prop*100
dta$year <- factor(dta$year)

m <- lm(lobby.pct ~ city.density + gov.partyR + house.rep.prop + divided.gov  + 
          house.np.diff + prop.state.transfers + year, data = dta)
summary(m)

N <- rep(c(dim(m$model)[1]), 4)

modelcoef = c(m$coefficients[3], m$coefficients[4], m$coefficients[5],
              m$coefficients[6])
modelse = c(coef(summary(m))[, 2][3], coef(summary(m))[, 2][4], 
            coef(summary(m))[, 2][5], coef(summary(m))[, 2][6])

ylo = modelcoef - qt(.975, N)*(modelse)
yhi = modelcoef + qt(.975, N)*(modelse)
names = c("Republican Governor", "Proportion House Republicans",
          "Divided Government", "House Polarization")
dfplot = data.frame(names, modelcoef, modelse, ylo, yhi)

(p <- ggplot(dfplot, aes(names, modelcoef, ymin=ylo, ymax=yhi)) +
    geom_pointrange(size = .25) + 
    coord_flip() + geom_hline(yintercept = 0, lty=2) + xlab("") + ylab("") +
    theme_bw(base_family = "serif", base_size = 13)  +
    theme(plot.background = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank()))



##########################################################################################################
## Figure 4.6 State Transfers and City Lobbying Across and Within States
##########################################################################################################

dta <- read.csv("state_data.csv")

dta$gov.partyR <- ifelse(dta$gov.party == "R", 1, 0)
dta$leg.controlR <- ifelse(dta$leg.control == "R", 1, 0)
dta$house.rep.pct <- dta$house.rep.prop*100
dta$year <- factor(dta$year)

m1 <- lm(lobby.pct ~ city.density + gov.partyR + house.rep.prop + divided.gov  + 
           house.np.diff + prop.state.transfers + state + year, data = dta)

p1 <- ggplot(dta, aes(x = prop.state.transfers, y = lobby.pct)) +
    coord_cartesian(ylim=c(5,30)) +
    stat_summary_bin(geom = "point") +
    geom_smooth(method = lm, formula = y ~ poly(x, 1), 
                color = "black", se = FALSE) +
    labs(y = "% Cities Lobbying", 
         x = "State Transfers as % of City Budgets")  +
    theme_bw(base_family = "serif", base_size = 11) +
    theme(plot.background = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank()) 


## Residualize state transfers and prop. cities lobbying by state and year

m <- lm(prop.state.transfers ~ state + year, data = dta)
m1 <- lm(lobby.pct ~ state + year, data = dta[!is.na(dta$prop.state.transfers), ])

plot.dta <- data.frame(x = resid(m),
                       y = resid(m1))


p2 <- ggplot(plot.dta, aes(x = x, y = y)) +
    scale_x_continuous(limits = c(-5, 5), 
                       breaks = c(-5, -2.5, 0, 2.5, 5)) +
    stat_summary_bin(geom = "point") +
    geom_smooth(method = lm, formula = y ~ poly(x, 1), 
                color = "black", se = FALSE) +
    coord_cartesian(ylim=c(-4,4)) +
    labs(y = expression(paste(Delta, " % Cities Lobbying")), 
         x = expression(paste(Delta," State Transfers as % of City Budgets")))  +
    theme_bw(base_family = "serif", base_size = 11) +
    theme(plot.background = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank())


(fig <- ggarrange(p1, p2, ncol = 2, nrow = 1))


##########################################################################################################
## Figure 4.7 Within-State Effect of Lobbying on State Transfers As a Proportion of City Budgets
##########################################################################################################

dta <- read.csv("state_data.csv")

dta$year <- factor(dta$year)

m <- lm(prop.state.transfers ~ lobby.pct2 + lobby.pct1 + lobby.pct + lobby.pct.lag1 + state + year, data = dta)

## Format data for figure
beta <- c(m$coefficients[2], m$coefficients[3], m$coefficients[4], m$coefficients[5])
se <- c(coef(summary(m))[,2][2], coef(summary(m))[,2][3], coef(summary(m))[,2][4], coef(summary(m))[,2][5])
figdata <- data.frame(beta = c(beta), se = c(se), x = c(-1, 0, 1, 2))

## Figure: Leads and Lags
(p <- ggplot(figdata, aes(x = x, y = beta)) + 
    geom_point(size = 2) +
    geom_errorbar(aes(ymin = beta - (1.6 * se), ymax = beta + (1.6 * se)), 
                  size = .5, width = .1) +
    geom_hline(yintercept = 0, linetype = 2) +
    scale_x_continuous(breaks = c(-1, 0, 1, 2), 
                       labels = c("Year t-2", "Year t-1", "Year t", "Year t+1")) +
    geom_vline(xintercept = .5) +
    theme_bw(base_family = "serif", base_size = 12)  +
    labs(y = "Effect on State Transfers as % of City Budgets", 
         x = "% Cities Lobbying") +
    theme(plot.background = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank()))


